Modelling E. coli die-off in mini trials

change log

• 2019-05-03: fixed sCOD values, fixed methane values in export
• 2019-05-02: after meeting, update calculations for day 0 values.
• 2019-05-01: setup .Rmd, incorporate inital day 0 values

Introduction

The Data

Main dataset

Data is stored in this repo’s data dir. The raw directory contains the raw data files recieved, which were copied to the clean dir, and modified as needed.

The bulk of the data was from file “Box Behnken analysis.xlsx”, recieved from RN via Skype on 2019-02-27, and this file was an updated version of the earlier “All Raw Data Mini Trial.xlsx” file. The data were in the “Block on rep” and “methane” tabs.

The main changes to “Box Behnken analysis” that were made: - label OLR was changed to OL to reflect the time-independant nature of the value - label HRT was changed to time, for the inverse reason - the methane aggretation was fixed to prevent values from being summed twice - addition of an “export” tab, for exporting to csv.

All these changes can be found in “./data/clean/Box Behnken analysis_nick_27.02.19.xlsx”, and exported as “./data/clean/trial_data.csv”.

Day 0 Values

After discussion about the experimental setup and the move to a time-based model rather than the Box Behnken analysis, we requested the Day 0 values for the 3 feedstocks, innoculum, and their mixtures. The file “Mini Trial FS FIB TS&VS.xlsx” was recieved 2019-04-30 from SN.

Data was aggregated and exported as “./data/clean/day0.csv” from the “day0” worksheet I added.

Processed data

The results of the two datasets are combined and written to the clean folder as cleaned_combined_minitrial_data.csv.

A note about experimental design and terms

Three things are mixed in these reractors: slurry, Food waste, and the iccoculum (from an AD reactor) Recipe may be used interchangably with Feedstock. There are 3 recipes, with different slurry:food waste ratios of 1:2, 1:3, and 3:1.

• OL: organic load, expressed in g VS/L
• sCOD: the soluble Chemical Oxygen Demand

Todo

• incorporate innoculumn sCOD values when recieved; either retest orr impute pH.

Data explorations

First thing we did was read in the main data file and merge it with the day 0 data. Because the day zero are true for all temperatures, we duplicate the data for the 19 and 55 degree conditions. Then, we convert log-scale the the indicator bacteria counts, adding a dummy value of 1 to prevent errors of logging 0.

main_data <- read.csv("./data/clean/trial_data.csv")
day0 <- read.csv("./data/clean/day0.csv", skip=3, skipNul = T, stringsAsFactors = F)
day0 <- day0 %>% select(-fs_conc, -fs_vol, -inoc_conc, -inoc_vol) %>% spread(key = measurement, value = Value)
day0<- rbind(
  day0,
  day0 %>% transform(Temp=replace(Temp,Temp==37, 19)),
  day0 %>% transform(Temp=replace(Temp,Temp==37, 55))
) 
names(main_data)

##  [1] "Temp"        "OL"          "Time"        "Recipe"      "Ammonia"    
##  [6] "pH"          "sCOD"        "TS"          "VS"          "Coliforms"  
## [11] "Ecoli"       "Enterococci" "Methane"

names(day0)

##  [1] "Recipe"      "Temp"        "OL"          "Time"        "Ammonia"    
##  [6] "Coliforms"   "Ecoli"       "Enterococci" "Methane"     "pH"         
## [11] "sCOD"        "TS"          "VS"

sn_raw <- rbind(main_data, day0)

sn <- sn_raw %>%
  rowwise() %>% 
  transform(
    Recipe=as.factor(Recipe),
    Ecoli =  log10(Ecoli + 1),
    Coliforms = log10(Coliforms + 1),
    Enterococci = log10(Enterococci + 1))
summary(sn)

##       Temp          OL           Time         Recipe      Ammonia    
##  Min.   :19   Min.   :0.5   Min.   :0.000   FWS12:69   Min.   : 865  
##  1st Qu.:19   1st Qu.:0.5   1st Qu.:3.000   FWS13:69   1st Qu.:1452  
##  Median :37   Median :2.0   Median :6.000   FWS31:69   Median :1625  
##  Mean   :37   Mean   :2.0   Mean   :5.217              Mean   :1619  
##  3rd Qu.:55   3rd Qu.:3.5   3rd Qu.:6.000              3rd Qu.:1768  
##  Max.   :55   Max.   :3.5   Max.   :9.000              Max.   :2305  
##        pH             sCOD             TS              VS       
##  Min.   :6.840   Min.   : 6522   Min.   :4.840   Min.   :3.120  
##  1st Qu.:7.240   1st Qu.: 8262   1st Qu.:5.570   1st Qu.:3.410  
##  Median :7.440   Median :10111   Median :5.740   Median :3.560  
##  Mean   :7.436   Mean   :10502   Mean   :5.874   Mean   :3.705  
##  3rd Qu.:7.510   3rd Qu.:12261   3rd Qu.:6.035   3rd Qu.:3.860  
##  Max.   :8.060   Max.   :22965   Max.   :8.095   Max.   :5.371  
##    Coliforms         Ecoli        Enterococci       Methane      
##  Min.   :0.000   Min.   :0.000   Min.   :0.000   Min.   :  0.00  
##  1st Qu.:2.249   1st Qu.:2.154   1st Qu.:4.511   1st Qu.:  5.75  
##  Median :4.772   Median :4.502   Median :5.178   Median : 19.50  
##  Mean   :3.857   Mean   :3.675   Mean   :4.757   Mean   : 40.89  
##  3rd Qu.:5.380   3rd Qu.:5.255   3rd Qu.:5.871   3rd Qu.: 75.55  
##  Max.   :8.020   Max.   :8.020   Max.   :7.322   Max.   :188.30

It may prove useful to note the components of the three recipes. The organic loading rates were calcualted based on the volitile solid quantity. the recepies are ratios of 2 of three ingrrerdiates (not including the innoculum, which is added as discussed above): dairy cattle slurry(DCS), Fats Oils and Grease (FOG), and restauraant food waste (RFW). - FWS12: 1:2 FOG:DCS - FWS13: DCS:RFW 1:3 - FWS31: DCS:RFW 3:1

We add those proportion to the data. Then, we save the cleaned, prepated data.

recipes <- data.frame(
  Recipe= c("FWS12","FWS13","FWS31"),
  RFW= as.factor(round(c(0, 3/4, 1/4), 2)),
  DCS=as.factor(round(c(1/3, 1/4, 3/4),2)),
  FOG=as.factor(round(c(2/3, 0, 0), 2))
)
sn <- merge(sn, recipes, by="Recipe")
write.csv(sn, "data/clean/cleaned_combined_minitrial_data.csv")

Note that the conditions were set to have 3 different loading rates, which we will treat as categorical because while they are numerical values of .5, 2, or 3.5, these are not observed values. There for we set them to be factors (which allows for categories to be ranked by a inheirent value).

Next, we can do some explortatory analysis:

sn$OL <- as.factor(sn$OL)
ggscatmat(sn,  color="Recipe", alpha = .3) + theme(legend.position = "bottom") 

## Warning in ggscatmat(sn, color = "Recipe", alpha = 0.3): Factor variables
## are omitted in plot

#ggpairs(sn, aes(color=Recipe, alpha=.5))

Initial observations: - the original Box Behnken design can be seen looking at the top left area, seeing the distribution of values between Time, OL, and Temp. - E. coli and Coliforms are, unsurprisingly, very highly correlated. - TS and VS are highly correlated

Lets replot without some of the highly correlated ones that

ggscatmat(sn %>% select(-TS, -Coliforms),  color="Recipe", alpha = .3) + theme(legend.position = "bottom")

## Warning in ggscatmat(sn %>% select(-TS, -Coliforms), color = "Recipe",
## alpha = 0.3): Factor variables are omitted in plot

Lets look at the indocators ### Enterococci

ggplot(sn, aes(x=Time, y=Enterococci, color=as.factor(Temp))) +
  geom_point() + facet_wrap(~OL+Recipe)  + geom_smooth(method = lm) + 
  labs(caption=fn("Effect of OL and Recipe"), color="")

It looks like for FWS12 the levels of enterococci stay the same, but in all the FSW13, FWS31, the levels seems to increase or decrease over time. We can try to determine if that is due to

ggplot(sn, aes(x=Time, y=Enterococci, color=as.factor(Temp))) + 
  geom_point() + 
  labs(caption=fn("Effect of OL vs DCS"), color="") +
  facet_wrap(~OL+DCS)  + 
  geom_smooth(method = lm)

ggplot(sn, aes(x=Time, y=Enterococci, color=as.factor(Temp))) + 
  labs(caption=fn("Effect of OL vs RFW"), color="") +
  geom_point() + facet_wrap(~OL+RFW)  + 
  geom_smooth(method = lm)

ggplot(sn, aes(x=Time, y=Enterococci, color=as.factor(Temp))) +
  geom_point() + 
  labs(caption=fn("Effect of Temp"), color="") +
  facet_wrap(~RFW)  + geom_smooth(method = lm)

ggplot(sn, aes(x=Temp, y=Enterococci, color=as.factor(RFW))) + 
  labs(caption=fn("Effect of RFW"), color="") +
  geom_point() + geom_smooth(method = lm)

Those rates seem to indicate 2 things:

• that the higher proportion of restaurant food waste, the better enterococci survive
• regardless of the proportion of RFW, increased temperature inhibits Enterococci

E. coli

ggplot(sn, aes(x=Time, y=Ecoli, color=as.factor(Temp))) +
  geom_point() + 
  facet_wrap(~OL+Recipe)  +
  geom_smooth(method = lm) + 
  labs(caption=fn("Effect of OL and Recipe"), color="")

E coli counts seem to depend largely on Temp rather than OL or Recipe.

ggplot(sn, aes(x=Time, y=Ecoli, color=as.factor(Temp))) +
  geom_point()  + geom_smooth(method = lm) + 
  labs(caption=fn("Effect of Temperature"), color="")

This may indicate that while die-off occurs at cold temperature, higher temperatue increases the die-off rate.

Coliforms

Coliforms should be have the same as E. coli (as we saw theirr values being highly correlated), but lets observe:

ggplot(sn, aes(x=Time, y=Coliforms, color=as.factor(Temp))) +
  geom_point()  + geom_smooth(method = lm) + 
  labs(caption=fn("Effect of Temperature"), color="")

Coliform counts too seem to depend largely on Temp rather than OL or Recipe.

Ammonia

ggplot(sn, aes(x=Time, y=Ammonia, color=as.factor(Temp))) +
  geom_point()  + geom_smooth(method = lm) + 
  facet_wrap(~Recipe )  +
  labs(caption=fn("Effect of Temperature and Recipe"), color="")

Across all recipes, Ammonia decreases as a function of time

Methane

Methane is the only thing leaving the bioreactors (as biogas), and therefore was measured daily until detructive sampling at 3, 6, or 9 days. We have aggregated the data in the original excel sheet.

ggplot(sn, aes(x=Time, y=Methane, color=as.factor(Recipe))) +
  geom_point()  + geom_smooth(method = lm) + 
  facet_wrap(~Temp )  +
  labs(caption=fn("Effect of Temperature and Recipe"), color="")

ggplot(sn, aes(x=Time, y=Methane, color=as.factor(Temp))) +
  geom_point()  + geom_smooth(method = lm) + 
  facet_wrap(~Recipe )  +
  labs(caption=fn("Effect of Temperature and Recipe"), color="")

ggplot(sn, aes(x=Time, y=Methane, color=as.factor(Temp))) +
  geom_point()  + geom_smooth(method = lm) + 
  facet_wrap(~Recipe + OL )  +
  labs(caption=fn("Effect of Temperature and OL"), color="")

It looks like the rate of poduction is highest at 37, and that at the optimal feed stocks with a high OL yield high prroduction, expecially at 37.

pH

ggplot(sn, aes(x=Time, y=pH, color=as.factor(Recipe))) +
  geom_point()  + geom_smooth(method = lm) + 
  facet_wrap(~Temp )  +
  labs(caption=fn("Effect of Temperature and Recipe"), color="")

The range of pH is rerlatively tight, but it appears to drop with time at the temperature extremes, staying the most stable at 37

Modelling

stan_data = list(
    nbugs=3,
    nrows=nrow(sn),
    Temp=(sn$Temp -19)/18 + 1,
    Ammonia=sn$Ammonia,
    pH=sn$pH,
    Ecoli=sn$Ecoli,
    Coliforms=sn$Coliforms,
    Time=sn$Time,
    OL=ifelse(sn$OL==.5, 1, ifelse(sn$OL==3.5, 3, sn$OL)),
    Enterococci=sn$Enterococci,
    sCOD=sn$sCOD,
    VS=sn$VS,
    TS=sn$TS,
    Recipe=as.numeric(sn$Recipe),
    Methane=sn$Methane
  )
fit <- stan(
  file = "./models/model_min.stan",
  verbose = T,
  chains=4,
  data = stan_data
)

## 
## TRANSLATING MODEL 'model_min' FROM Stan CODE TO C++ CODE NOW.
## successful in parsing the Stan model 'model_min'.
## 
## CHECKING DATA AND PREPROCESSING FOR MODEL 'model_min' NOW.
## 
## COMPILING MODEL 'model_min' NOW.
## 
## STARTING SAMPLER FOR MODEL 'model_min' NOW.

## Warning: There were 1346 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

## Warning: Examine the pairs() plot to diagnose sampling problems

plot_pars <- names(fit)[!startsWith(prefix = "bugs", names(fit))]
plot(fit, pars=plot_pars, plotfun = "trace", inc_warmup = F)

mcmc_rhat(rhat(fit)[!startsWith(prefix = "bugs", names(fit))]) + yaxis_text(hjust = 0)

ggplot(data.frame(x = c(-5, 5)), aes(x)) +
    stat_function(fun = dcauchy, n = 1e3, args = list(location = 1, scale = 1))

stan_dens(fit, pars=c("recipe_effect"))

plot(fit, pars=c("recipe_effect"))

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

plot(fit, pars=c("temp_effect"))

## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

plot(fit, pars=c("lp__"))

## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)

#stan_model(fit)
(mle = optimizing(stan_model("./models/model_min.stan"), data=stan_data))

## $par
## recipe_effect[1] recipe_effect[2] recipe_effect[3]   temp_effect[1] 
##     1.786915e+00    -4.330820e-01    -6.778961e-01     4.311390e-02 
##   temp_effect[2]   temp_effect[3]     ol_effect[1]     ol_effect[2] 
##    -9.153030e-01     1.673147e-01     7.821730e+00     1.626072e+01 
##     ol_effect[3]            sigma      bugs_hat[1]      bugs_hat[2] 
##     1.787536e+01     2.160753e+30     1.183533e+05     1.125921e+05 
##      bugs_hat[3]      bugs_hat[4]      bugs_hat[5]      bugs_hat[6] 
##     1.657377e+05     1.254627e+05     7.855682e+04     1.099506e+05 
##      bugs_hat[7]      bugs_hat[8]      bugs_hat[9]     bugs_hat[10] 
##     1.748567e+05     1.679054e+05     2.132778e+05     1.667928e+05 
##     bugs_hat[11]     bugs_hat[12]     bugs_hat[13]     bugs_hat[14] 
##     1.571769e+05     2.121249e+05     1.706268e+05     1.724419e+05 
##     bugs_hat[15]     bugs_hat[16]     bugs_hat[17]     bugs_hat[18] 
##     6.044487e+04     6.027080e+04     1.651412e+05     1.685152e+05 
##     bugs_hat[19]     bugs_hat[20]     bugs_hat[21]     bugs_hat[22] 
##     2.209633e+05     8.191879e+04     7.832521e+04     1.726854e+05 
##     bugs_hat[23]     bugs_hat[24]     bugs_hat[25]     bugs_hat[26] 
##     1.166730e+05     1.688998e+05     1.087502e+05     1.694282e+05 
##     bugs_hat[27]     bugs_hat[28]     bugs_hat[29]     bugs_hat[30] 
##     7.786855e+04     4.105228e+05     1.130727e+05     1.178734e+05 
##     bugs_hat[31]     bugs_hat[32]     bugs_hat[33]     bugs_hat[34] 
##     1.116322e+05     1.810452e+05     8.631004e+04     1.725937e+05 
##     bugs_hat[35]     bugs_hat[36]     bugs_hat[37]     bugs_hat[38] 
##     5.615964e+04     1.154729e+05     1.537116e+05     5.951931e+04 
##     bugs_hat[39]     bugs_hat[40]     bugs_hat[41]     bugs_hat[42] 
##     1.860855e+05     1.237718e+05     1.786968e+05     1.685509e+05 
##     bugs_hat[43]     bugs_hat[44]     bugs_hat[45]     bugs_hat[46] 
##     1.690371e+05     5.822546e+04     1.644123e+05     1.365916e+05 
##     bugs_hat[47]     bugs_hat[48]     bugs_hat[49]     bugs_hat[50] 
##     2.109678e+05     1.952756e+05     1.583445e+05     1.165913e+05 
##     bugs_hat[51]     bugs_hat[52]     bugs_hat[53]     bugs_hat[54] 
##     8.122379e+04     1.608862e+05     1.305319e+05     6.703936e+04 
##     bugs_hat[55]     bugs_hat[56]     bugs_hat[57]     bugs_hat[58] 
##     1.332220e+05     8.295257e+04     6.266193e+04     2.928694e+05 
##     bugs_hat[59]     bugs_hat[60]     bugs_hat[61]     bugs_hat[62] 
##     2.928694e+05     4.105228e+05     7.382467e+04     1.021236e+05 
##     bugs_hat[63]     bugs_hat[64]     bugs_hat[65]     bugs_hat[66] 
##     2.928694e+05     4.105228e+05     6.844350e+04     1.021236e+05 
##     bugs_hat[67]     bugs_hat[68]     bugs_hat[69]     bugs_hat[70] 
##     1.748564e+05     1.366878e+05     1.021236e+05     1.419258e+05 
##     bugs_hat[71]     bugs_hat[72]     bugs_hat[73]     bugs_hat[74] 
##     1.123880e+05     2.021011e+05     1.378997e+05     7.966518e+04 
##     bugs_hat[75]     bugs_hat[76]     bugs_hat[77]     bugs_hat[78] 
##     8.177432e+04     1.809906e+05     1.855631e+05     8.421338e+04 
##     bugs_hat[79]     bugs_hat[80]     bugs_hat[81]     bugs_hat[82] 
##     2.449770e+05     2.028970e+05     2.518740e+05     1.412786e+05 
##     bugs_hat[83]     bugs_hat[84]     bugs_hat[85]     bugs_hat[86] 
##     2.450736e+05     7.424722e+04     2.154774e+05     7.825904e+04 
##     bugs_hat[87]     bugs_hat[88]     bugs_hat[89]     bugs_hat[90] 
##     1.328324e+05     3.188351e+05     8.410832e+04     1.836829e+05 
##     bugs_hat[91]     bugs_hat[92]     bugs_hat[93]     bugs_hat[94] 
##     1.436913e+05     3.188351e+05     5.950899e+04     1.406585e+05 
##     bugs_hat[95]     bugs_hat[96]     bugs_hat[97]     bugs_hat[98] 
##     2.409977e+05     8.668223e+04     1.860511e+05     1.419259e+05 
##     bugs_hat[99]    bugs_hat[100]    bugs_hat[101]    bugs_hat[102] 
##     7.309299e+04     9.629382e+04     1.944402e+05     2.048183e+05 
##    bugs_hat[103]    bugs_hat[104]    bugs_hat[105]    bugs_hat[106] 
##     1.746205e+05     1.679254e+05     2.450736e+05     1.309019e+05 
##    bugs_hat[107]    bugs_hat[108]    bugs_hat[109]    bugs_hat[110] 
##     2.362202e+05     5.950849e+04     2.317106e+05     1.465055e+05 
##    bugs_hat[111]    bugs_hat[112]    bugs_hat[113]    bugs_hat[114] 
##     2.159638e+05     1.877536e+05     1.717681e+05     1.474174e+05 
##    bugs_hat[115]    bugs_hat[116]    bugs_hat[117]    bugs_hat[118] 
##     7.355702e+04     2.017443e+05     1.870239e+05     2.498922e+05 
##    bugs_hat[119]    bugs_hat[120]    bugs_hat[121]    bugs_hat[122] 
##     2.351592e+05     2.450736e+05     5.840009e+04     9.629382e+04 
##    bugs_hat[123]    bugs_hat[124]    bugs_hat[125]    bugs_hat[126] 
##     1.523874e+05     6.468478e+04     1.262270e+05     7.713658e+04 
##    bugs_hat[127]    bugs_hat[128]    bugs_hat[129]    bugs_hat[130] 
##     9.629382e+04     1.391061e+05     2.322411e+05     1.756472e+05 
##    bugs_hat[131]    bugs_hat[132]    bugs_hat[133]    bugs_hat[134] 
##     1.304193e+05     1.422432e+05     1.339400e+05     1.400714e+05 
##    bugs_hat[135]    bugs_hat[136]    bugs_hat[137]    bugs_hat[138] 
##     7.228446e+04     1.357281e+05     7.471261e+04     3.188351e+05 
##    bugs_hat[139]    bugs_hat[140]    bugs_hat[141]    bugs_hat[142] 
##     1.551771e+05     2.276386e+05     2.005904e+05     8.619428e+04 
##    bugs_hat[143]    bugs_hat[144]    bugs_hat[145]    bugs_hat[146] 
##     1.370592e+05     1.279139e+05     1.946281e+05     1.376996e+05 
##    bugs_hat[147]    bugs_hat[148]    bugs_hat[149]    bugs_hat[150] 
##     6.227995e+04     1.327276e+05     1.329684e+05     2.005906e+05 
##    bugs_hat[151]    bugs_hat[152]    bugs_hat[153]    bugs_hat[154] 
##     1.848310e+05     7.077057e+04     9.186715e+04     2.086660e+05 
##    bugs_hat[155]    bugs_hat[156]    bugs_hat[157]    bugs_hat[158] 
##     2.327458e+05     1.322461e+05     6.098742e+04     8.526667e+04 
##    bugs_hat[159]    bugs_hat[160]    bugs_hat[161]    bugs_hat[162] 
##     1.293572e+05     8.202186e+04     6.737795e+04     1.366524e+05 
##    bugs_hat[163]    bugs_hat[164]    bugs_hat[165]    bugs_hat[166] 
##     1.558638e+05     2.013722e+05     1.332089e+05     1.832947e+05 
##    bugs_hat[167]    bugs_hat[168]    bugs_hat[169]    bugs_hat[170] 
##     1.353303e+05     2.008834e+05     5.562575e+04     8.306325e+04 
##    bugs_hat[171]    bugs_hat[172]    bugs_hat[173]    bugs_hat[174] 
##     1.320042e+05     1.844476e+05     1.456484e+05     2.096003e+05 
##    bugs_hat[175]    bugs_hat[176]    bugs_hat[177]    bugs_hat[178] 
##     1.475041e+05     2.518917e+05     1.416325e+05     7.065319e+04 
##    bugs_hat[179]    bugs_hat[180]    bugs_hat[181]    bugs_hat[182] 
##     2.048193e+05     1.559704e+05     2.459603e+05     1.702286e+05 
##    bugs_hat[183]    bugs_hat[184]    bugs_hat[185]    bugs_hat[186] 
##     6.024708e+04     7.848485e+04     1.755901e+05     2.132718e+05 
##    bugs_hat[187]    bugs_hat[188]    bugs_hat[189]    bugs_hat[190] 
##     9.079141e+04     1.993639e+05     2.307861e+05     2.133430e+05 
##    bugs_hat[191]    bugs_hat[192]    bugs_hat[193]    bugs_hat[194] 
##     1.498855e+05     1.360102e+05     8.133175e+04     1.571327e+05 
##    bugs_hat[195]    bugs_hat[196]    bugs_hat[197]    bugs_hat[198] 
##     7.006891e+04     9.079141e+04     1.993639e+05     2.307861e+05 
##    bugs_hat[199]    bugs_hat[200]    bugs_hat[201]    bugs_hat[202] 
##     1.906687e+05     8.063218e+04     2.239904e+05     1.323661e+05 
##    bugs_hat[203]    bugs_hat[204]    bugs_hat[205]    bugs_hat[206] 
##     1.387446e+05     1.433174e+05     9.079141e+04     1.993639e+05 
##    bugs_hat[207] 
##     2.307861e+05 
## 
## $value
## [1] -14603.89
## 
## $return_code
## [1] 0